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Abstract 

We consider simulation of spatially one-dimensional space-time 
fractional diffusion. Whereas in an earlier paper of ours (Eur. 
Phys. J. Special Topics, Vol. 193, 119-132 (2011); E-print: 



http://arxiv.org/abs/1104.4041j, we have developed the basic 



theory of what we call parametric subordination via three-fold split- 
ting applied to continuous time random walk with subsequent passage 
to the diffusion limit, here we go the opposite way. Via Fourier-Laplace 
manipulations of the relevant fractional partial differential equation of 
evolution we obtain the subordination integral formula that teaches 
us how a particle path can be constructed by first generating the 
operational time from the physical time and then generating in 
operational time the spatial path. By inverting the generation of 
operational time from physical time we arrive at the method of 
parametric subordination. Due to the infinite divisibility of the stable 
subordinator we can simulate particle paths by discretization where 
the generated points of a path are precise snapshots of a true path. 
By refining the discretization more and more fine details of a path 
become visible. 



^E-print based on an updated version of Chapter 10, pp 227-261 in S.C. Lim, J. 
Klafter and R. Metzler (Editors): Fractional Dynamics, Recent Advances, World Scientific, 
Singapore 2012 (ISBN 9814340588, 9789814340588), 
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1 Introduction 



The purpose of this chapter is to describe our method of parametric 
subordination to produce particle trajectories for the so-called fractional 
diffusion processes. 

By replacing in the common diffusion equation the first order time derivative 
and the second order space derivative by appropriate fractional derivatives we 
obtain a fractional diffusion equation whose solution describes the temporal 
evolution of the density of an extensive quantity, e.g. of the sojourn 
probability of a diffusing particle. 

After giving a survey on analytic methods for determination of the solution 
(this is the macroscopic aspect) we turn attention to the problem of 
simulation of particle trajectories (the microscopic aspect). By some authors 
such simulation is called "particle tracking", see e.g. |65j. 

As an approximate method among physicists the so-called Continuous Time 
Random Walk {CTRW) is very popular. On the other hand, it is possible 
to produce a sequence of precise snapshots of a true trajectory. This is 
achieved by a change from the "physical time" to an "operational time" in 
which the simulation is carried out. By two Markov processes happening in 
operational time the running of physical time and the motion in space are 
produced. Then, elimination of the operational time yields a picture of the 
desired trajectory. It is remarkable that so by combination of two Markov 
processes a non-Markovian process is generated. 

The two Markov processes can be obtained and analyzed in two ways: 

(a) from the CTRW model by a well-scaled passage to the " diffusion limit" , 

(b) directly from an integral representation of the fundamental solution of 
the fractional diffusion equation. 

We have developed way (a) in our 2007 paper [IB] via passage to the diffusion 
limit in the Cox- Weiss solution formula for CTRW and by the technique 
of splitting the CTRW into three separate walks and passing in each 
of these to the diffusion limit in our more recent recent paper [17]. For 
another access (more oriented towards measure-theoretic theory of stochastic 
processes) see the recent papers [TSl [SO]- In [20] the authors also treat 
the problem of subordination for diffusion with distributed orders of time- 
fractional differentiation. 
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The plan of our chapter is as follows. In Section 2 we provide for the reader's 
convenience some preliminary notions and notations as a mathematical 
background for our further analysis. In Section 3 we introduce the space- 
time fractional diffusion equation, based on the Riesz-Feller and Caputo 
fractional derivatives, and we present the fundamental solution.. In Section 4 
we provide the stochastic interpretation of the space-time fractional diffusion 
equation discussing the concepts of subordination, the main goal of this 
chapter. Finally, in Section 5 we show some graphical representations along 
with conclusions. 



2 Notions and Notations 

In this Section we survey some preliminary notions including Fourier and 
Laplace transforms, special functions of Mittag-LefHer and Wright type and 
Levy stable probability distributions. 

Since in what follows we shall meet only real or complex-valued functions 
of a real variable that are defined and continuous in a given open interval 
I = (a, b) , — oo < a < b < +oo , except, possibly, at isolated points where 
these functions can be infinite, we restrict our presentation of the integral 
transforms to the class of functions for which the Riemann improper integral 
on / absolutely converges. In so doing we follow Marichev [33j and we denote 
this class by L'^{I) or L'^{a, b) . 



2.1 The Fourier transform 

Let 

/ + 00 
e+^^^ f{x)dx, i^eR, (2.1a) 
-oo 

be the Fourier transform of a function f{x) G L'^(IR), and let 

f{x)=T-'[f{Ky,xj = — J e-'^^f{K)dK, xem, (2.16) 

be the inverse Fourier transform^. 

■^If f{x) is piecewise differentiable, then the formula (2.1b) holds true at all points 
where f{x) is continuous and the integral in it must be understood in the sense of the 
Cauchy principal value. 
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Related to the Fourier transform is the notion of pseudo-differential operator. 
Let us recall that a generic pseudo-differential operator A, acting with respect 
to the variable a; G IR , is defined through its Fourier representation, namely 

/ + 00 
e'^^A[fix)]dx = AiK)fiK), (2.2) 
-oo 

where A[k) is referred to as symbol of A , formally given as 
A{k) = (Ae"*'^^) e+'^\ 

2.2 The Laplace transform 

Let 

^ /"OO 

f{s) = C{f{ty,s}= e-^^f{t)dt, Re(s)>a;, (2.3a) 
Jo 

be the Laplace transform of a function f{t) G £^(0,T) , VT > and let 

f{t) = C-Hf{s)-t\ = — e'^f{s)ds, Re(s)=7>a/, (2.36) 

I J 2m J^_i^ 

with t > , be the inverse Laplace transforn|^ 

2.3 The auxihary functions of Mittag-Lefiier type 

The Mittag-Leffler functions, that we denote by Ea{z), Ea^p{z) are so named 
in honour of Gosta Mittag-Leffler, the eminent Swedish mathematician, who 
introduced and investigated these functions in a series of notes starting from 
1903 in the framework of the theory of entire functions[lSl HSl 113 • The 
functions are defined by the series representations, convergent in the whole 
complex plane (D 

n=0 

■^A sufficient condition of tlie existence of tlie Laplace transform is that tlie original 
function is of exponential order as i — > oo . This means that some constant a / exists such 
that the product e~°^* 1/(^)1 is bounded for all t greater than some T . Then f{s) exists 
and is analytic in the half plane Re(s) > a/ . If /(i) is piecewise differentiable, then the 
formula (2.3b) holds true at all points where f{t) is continuous and the (complex) integral 
in it must be understood in the sense of the Cauchy principal value. 
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E^A^):=J2 \n. , Re(a)>0, /3gC. (2.5) 

Originally Mittag-Leffler assumed only the parameter a and assumed it as 
positive, but soon later the generalization with two complex parameters was 
considered by Wiman[60]. In both cases the Mittag-Leffler functions are 
entire of order l/Re(Q;). Generally Ea^i{z) = Ea{z). 

Using their series representations it is easy to recognize 

E,,(z) = EAz)=e', E,o(z) = ~ ^ 



E2,i{z^) = cosh(z) , E2,i{-z^) = los{z) , (2.6) 

i^..(.^) = ^, E2A-z') = '^. 



and more generally 

E^,f,{z) + E^^^{-z) =2E2aAz^) 

Ea^z) - Ea^-Z) = 2zE2a,a+0{z 



2^ 



(2.7) 



We note that in Chapter 18 of Vol. 3 of the handbook of the 1955 Bateman 
Project [5] devoted to Miscellaneous Functions, we find a valuable survey of 
these functions, which later were recognized as belonging to the more general 
class of Fox if-functions introduced after 1960. 

For our purposes relevant roles are played by the following auxiliary functions 
of the Mittag-Leffler type on support IR"*" defined as follows, where A > 0, 
along with their Laplace transforms 

e„(t;A):=K(-An-^^, (2.8) 

S -r A 

CaAt; A) := t^-' E^,, (-A n - 4^ , (2.9) 

S + A 

e„,„(t; A) := r-'E^^^ (-At°) = ^e,(-At") = - - , (2.10) 

Here we have used the sign for the juxtaposition of a function depending 
on t with its Laplace transform depending on s. Later we use this sign also 
for juxtaposition of a function depending on x with its Fourier transform 
depending on k. 
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Remark: We outline that the above auxihary functions (for restricted values 
of the parameters) turn out to be completely monotone (CM) functions so 
that they enter in some types of relaxation phenomena of physical relevance. 
We recall that a function f{t) is CM in IR + if (-l)"/"(t) > 0. The 
function e~* is the prototype of a CM function. For a Bernstein theorem, 
more generally they are expressed in terms of a (generalized) real Laplace 
transform of a positive measure 

POO 

f{t)= e-'' K{r)dr , K{r)>0. (2.11) 
Jo 

Restricting attention to the auxiliary function in two parameters, we can 
prove for A > that 

e^^f^{t;X) ■.= t^-^ Ea^p{-Xt") CM iff < a < /3 < 1 . (2.12) 

Using the Laplace transform we can prove, following Gorenflo and Mainardi 
[13], that for < a < 1 and A = 1 

1 t 0+ 

EA-n-^{ ",£F+^^,3F+^''' ' (2.13) 

r(l-a)"r(l-2a)'" t^+^^ 
and 

POO 

E„(-t°)= / e-'' K^{r) dr (2.14) 
Jo 

with 

K.{r) = i , r'^-'Ma-) ^ l_ ^ _ (2 ,5) 

TT r^" + 2 r" cos(a;7r) + 1 vr r r° + 2 cos(a7r) + r " 

2.4 The auxiliary functions of the Wright type 

The Wright function, that we denote by W^A.^tX-^)? is so named in honour of 
E. Maitland Wright, the eminent British mathematician, who introduced 
and investigated this function in a series of notes starting from 1933 in 
the framework of the asymptotic theory of partitions, see [611 [62| 163] . The 
function is defined by the series representation, convergent in the whole z- 
complex plane C, 

oo n 

W^,,{z):=Y, M W A>-l,/iGC. (2.16) 

nl r An + u 

n=0 ^ ' 
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Originally, Wright assumed A > 0, and, only in 1940 [M], he considered 
— 1 < A < 0. We note that in Chapter 18 of Vol. 3 of the handbook of the 
1955 Bateman Project |5j devoted to Miscellaneous Functions, we find an 
earlier analysis of these functions, which, similarly with the Mittag-Leffler 
functions, were later recognized as belonging to the more general class of Fox 
iif-functions introduced after 1960. However, in that Chapter, presumably 
for a misprint, the parameter A of the Wright function is restricted to be non 
negative. It is possible to prove that the Wright function is entire of order 
1/(1 + A) , hence it is of exponential type only if A > 0. For this reason 
we propose to distinguish the Wright functions in two kinds according to 
A > {first kind) and — 1 < A < {second kind). Both kinds of functions are 
related to the Mittag-Leffler function via Laplace transform pairs: in fact we 
have, see for details the appendix F of the recent book by Mainardi [28], for 
the case A > (Wright functions of the first kind) 

Wx,,{±r) - - ^A,M (±-] , A > , (2.17) 



s \ s ^ 

and for the case A = —u E (—1, 0) (Wright functions of the second kind), 

iy_,,^(-r) - S.,^+.(-s) , < z/ < 1 . (2.18) 

For our purposes relevant roles are played by the following auxiliary functions 
of the Wright type (of the second kind) 

F,{z) := W.,,^{-z) = V . , < z/ < 1 , (2.19) 

71=1 ^ ' 

and 

M,{z) := iy_.,i-.(-z) = ^ , < z/ < 1 , (2.20) 

f^^n\r[-un + {l-u)] 

interrelated through F,y{z) = v z My{z). The relevance of these functions was 
pointed out by Mainardi in his former analysis of the time fractional diffusion 
equation via Laplace transform. Restricting our attention to the M- Wright 
functions on support IR"^ we point out the Laplace transforms pairs 

M^{r) ^ E^{-s) , 0<z/<l. (2.21) 
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0<z/<l. (2.22) 

< z/ < 1 . (2.23) 

It was also proved in ^28j that the M- Wright function on support IR"*" is a 
probabihty density function (pdf) )that in the hterature is sometimes known 
as the density related to Mittag-Leffler probability distribution. Its absolute 
moments of order 6 > —1 in IR"*" are finite and turn out to be 

/ r'M,{r)dx = -^—^, 6>-l, 0<u<l. (2.24) 

Jo r(z/(^ + 1) 

We point out that in the limit u ^ 1^ the function Mi,{r), for r G IR"*", tends 
to the Dirac generalized function 6{r — 1). 

For our next purposes it is worthwhile to introduce the function in two 
variables 

M^{x,t) ■.= t-'' M^ixt-") , 0<z/<l, x,tGlR+, (2.25) 

which defines a spatial probability density in x evolving in time t with self- 
similarity exponent H = v. Of course for s G IR we can consider the 
symmetric version obtained from (2.25) multiplying by 1/2 and replacing x 
by Hereafter we provide a list of the main properties of this density, which 
can be derived from the Laplace and Fourier transforms for the corresponding 
Wright M-function in one variable. 

From Eq. (2.23) we derive the Laplace transform of My{x,t) with respect to 
t G IR + , 

C{M^{x,t)]t-^ s} = s''~^e~^^\ (2.26) 

From Eq. (2.21) we derive the Laplace transform of My{x,t) with respect to 
X G 1R + , 

C {M^(x, t)]x^ s} = E^ {-St") . (2.27) 

From the recent book by Mainardi[28] we recall the Fourier transform of 
M^{\x\,t) with respect to a; G IR, 

{M,{\x\,ty,x ^ k} = {-kH^") , (2.28) 



(1/r^) 



i-M.(l/0 - - 



—s 
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and, in particular, 



cos(kx) Mj,(x, t) dx = E2u,i{—i^'^ t"^") 
sin(Kx) M^{x, t) dx = nt" £'2^,1+1/ (—/t^ t'^"'^ 



It is worthwhile to note that for z/ = 1/2 we recover the Gaussian density 
evolving with time with variance = 2t 

2.5 The Levy stable distributions 

The term stable has been assigned by the French mathematician Paul Levy, 
who, in the twenties of the last century, started a systematic research in 
order to generalize the celebrated Central Limit Theorem to probability 
distributions with infinite variance. For stable distributions we can assume 
the following Definition: // two independent real random variables with 
the same shape or type of distribution are combined linearly with positive 
coefficients and the distribution of the resulting random variable has the same 
shape, the common distribution (or its type, more precisely) is said to be 
stable. 

The restrictive condition of stability enabled Levy (and then other authors) 
to derive the canonic form for the characteristic function of the densities of 
these distributions. Here we follow the parameterization by Feller [6] revisited 
in |T4] and in [29j. Denoting by i^^(^) ^ generic stable density in IR, where 
a is the index of stability and and 6 the asymmetry parameter, improperly 
called skewness, its characteristic function reads: 

Liix) - LliK) = exp [-<(«:)] , = |«:|« e^(sign «:)^vr/2 

< a < 2, 1^1 < min{a,2 - a} . 

We note that the allowed region for the real parameters a and 6 turns 
out to be a diamond in the plane {a, 9} with vertices in the points 
(0,0) , (1, 1) , (1, —1) , (2,0), that we call the Feller- Takayasu diamond, see 
Figure [1] For values of 6 on the border of the diamond (that is ^ = ±a if 
< a < 1, and 6 = ±(2 — a) if 1 < a < 2) we obtain the so-called extremal 
stable densities. 
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Figure 1: The Feller- Takayasu diamond. 



We note the symmetry relation x) = L^^(x), so that a stable density 
with 6* = is symmetric. 

Stable distributions have noteworthy properties on which the interested 
reader can be informed from the relevant existing literature. Here-after we 
recall some peculiar Properties: 

- Each stable density possesses a domain of attraction, see e.g. [6]. 

- Any stable density is unimodal and indeed bell-shaped, i.e. its n-th. derivative 
has exactly n zeros in H, see [TO] . 

- The stable distributions are self-similar and infinitely divisible. 

These properties derive from the canonic form (2.31) through the scaling 
property of the Fourier transform. 
Self- similarity means 

t) - exp ^ t) = t-'l^ L^ixlt'l^)] , (2.32) 

where t is a positive parameter. If t is time, then L^^{x, t) is a spatial density 
evolving in time with self-similarity. 

Infinite divisibility means that for every positive integer n, the characteristic 
function can be expressed as the nth power of some characteristic function, 
so that any stable distribution can be expressed as the n-fold convolution 
of a stable distribution of the same type. Indeed, taking in (2.31) 9 = 0, 
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without loss of generality, we have 

e-l'^l" = ^ Ll{x,t) = , (2.33) 

where 

[LO (x, t/n)] ™ := L° (x, t/n) * L° (x, t/n) * ■ ■ ■ * (x, t/n) (2.34) 

is the multiple Fourier convolution in IR with n identical terms. 

Only in special cases the inversion of the Fourier transform in (2.31) can 
be carried out using standard tables, and provides well-known probability 
distributions. 

For a = 2 (so ^ = 0), we recover the Gaussian pdf, that turns out to be 
the only stable density with finite variance, and more generally with finite 
moments of any order 6 > 0. In fact 

L°(x) = -l=e-^V4. (2.35) 
2a/ vr 

All the other stable densities have finite absolute moments of order S G 
[—1,0;) as we will later show. 

For a = 1 and |^| < 1, we get 

Te(r) = - r2 3fi) 

' TT [x + sin(07r/2)]2 + [cos(^7r/2)]2 ' ^ ' 

which for 6 = includes the Cauchy-Lorentz pdf, 



L?(x) = --^. (2.37) 

TT 1 + X"^ 



In the limiting cases 6 = ±1 for a = 1 we obtain the singular Dirac pdf's 

Lf\x) =6{x±l) . (2.38) 



In general, we must recall the power series expansions provided in [6]. We 
restrict our attention to x > since the evaluations for x < can be obtained 
using the symmetry relation. The convergent expansions of L^^ix) (x > 0) 
turn out to be; 
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for < a < 1 , |6'| < a : 

n=l 

for 1 < a < 2, \e\<2-a : 



(^) = — - ") • (2-40) 

n=l 



From the series in (2.39) and the symmetry relation we note that the extremal 
stable densities for < a < 1 are unilateral, precisely vanishing for x > if 
6 = a, vanishing for x<Oif^ = — a. In particular the unilateral extremal 
densities L~"(x) with < a < 1 have support H"'" and Laplace transform 
exp (— s"). For a = 1/2 we obtain the so-called Levy-Smirnov pdf: 

—3/2 

L-f(x) = |-^e-V(4x)^ x>0. (2.41) 

It is worth to note that the Gaussian pdf (2.35) and the Levy-Smirnov pdf 
(2.41) are well known in the treatment of the Brownian motion: the former 
as the spatial density on an infinite real line, the latter as the first passage 
time density on a semi-infinite line, see e.g. [6]. 

As a consequence of the convergence of the series in (2.39)-(2.40) and of the 
symmetry relation we recognize that the stable pdf^s with 1 < a < 2 are 
entire functions, whereas with < a < 1 have the form 

, , , [l/x] <l>i(x-°) for X > 0, 
(x) = C ' ' ^ ' 2.42 
^ ^ ^ (l/|x|)$2(|a;|-") forx < 0, ^ ^ 

where ^i{z) and $2(2;) are distinct entire functions. The case a = 1 with 
\9\ < 1 must be considered in the limit for a — )■ 1 of (2.39)-(2.40), because the 
corresponding series reduce to power series akin with geometric series in 1 /x 
and X, respectively, with a finite radius of convergence. The corresponding 
stable densities are no longer represented by entire functions, as can be noted 
directly from their exphcit expressions (2.36)-(2.37). 

From a comparison between the series expansions in (2.39)-(2.40) and in 
(2.19)-(2.20), we recognize that for x > our auxiliary functions of the 
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Wright type are related to the extremal stable densities as follows, see|32]. 

L-"(x) = -F„(a;-") = ^M,(x-"), 0<a<l, (2.43) 

L'^^-^ix) = - = - Mi/,(x) , 1< « < 2 . (2.44) 

a; a 

In the above equations, for a = 1, the skewness parameter turns out to be 
= — 1, so we get the singular hmit 

Ly\x) = Mi{x) = 6{x - 1) . (2.45) 

We do not provide here the asymptotic representations of the stable densities 
referring the interested reader to[29j. However, based on asymptotic 
representations, we can state the following: For < a < 2 the stable densities 
exhibit fat tails in such a way that their absolute moment of order 6 is finite 
only if — 1 < 5 < a. More precisely, one can show that for non-Gaussian, not 
extremal, stable densities the asymptotic decay of the tails is 

L^(x) = O , x^±oo. (2.46) 

For the extremal densities with a ^ \ this is valid only for one tail (as 
\x\ — )■ oo\ the other (as |x| — )■ oo) being of exponential order. For 1 < a < 
2 the extremal prf/'s are two-sided and exhibit an exponential left tail (as 
X —7- — oo) if = +(2 — a) , or an exponential right tail (as x — t- +oo) if 
6* = — (2 — a) . Consequently, the Gaussian pdf is the unique stable density 
with finite variance. Furthermore, when < a < 1, the first absolute moment 
is infinite so we should use the median instead of the non-existent expected 
value in order to characterize the corresponding pdf . 

Let us also recall a relevant identity between stable densities with index a 
and 1/a (a sort of reciprocity relation) pointed out in|6], that is, assuming 

x > 0, 

^L?/Jx-")=Lr(x), l/2<a<l, r = a(0 + !)-!. (2.47) 

The condition 1/2 < a < 1 implies 1<1/q;<2. A check shows that Q* falls 
within the prescribed range |^*| < a if |0| < 2 — 1/a. We leave as an exercise 
for the interested reader the verification of this reciprocity relation in the 
limiting cases a = 1/2 and a = 1. For more details on Levy stable densities 
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we refer the reader to specialized treatises, as[Sl 1251 ESI ESI ESI EE], where 
different notations are adopted. We hke to refer also to the 1986 paper by 
Schneider [55], where he first provided the Fox if-function representation of 
the stable distributions (with a ^ 1) and to the 1990 book by Takayasu[56]. 
where he first gave the diamond representation in the plane {a, 6}. 

3 The Space-Time Fractional Diffusion 

We now consider the Cauchy problem for the (spatially one- dimensional) 
space-time fractional diffusion (STFD) equation. 

tD^%{x,t) = ^D^u{x,t) , u{x,0) = 6{x), xeJR, t>0, (3.1) 

where {a , 6 , (3} are real parameters restricted to the ranges 



Here denotes the Riesz-Feller fractional derivative of order a and 

skewness 6, acting on the space variable x, and t-Df denotes the Caputo 
fractional derivative of order /3, acting on the time variable t. We recall 
the definitions of these fractional derivatives based on their representation in 
the Fourier and Laplace transform domain, respectively. So doing we avoid 
the subtleties lying in the inversion of the corresponding fractional integrals, 
see e.g. the 2001 survey by Mainardi et al.[29j. For general information on 
fractional integrals and derivatives we recommend the books [251 SSI EI]- 

3.1 The Riesz-Feller space-fractional derivative 

We define the Riesz-Feller derivative as the pseudo-differential operator 
whose symbol is the logarithm of the characteristic function of a general Levy 
strictly stable probability density with index of stability a and asymmetry 
parameter 9 (improperly called skewness). As a consequence of Eq. (2.31), 
for a sufficiently well-behaved function /(x), we define the Riesz-Feller space- 
fractional derivative of order a and skewness 6 via the Fourier transform 



0<a<2, 1^1 < min{a,2 - a} , < (3 < 1 . 



(3.2) 




(3.3) 
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Notice that ^^sign/t _ j-^ (gign/t) 6-n/2]. For = we have a symmetric 
operator with respect to x , which can be interpreted as 



a/2 



as can be formally deduced by writing — = — (k^)"/^ . We thus recognize 
that the operator Dq is related to a power of the positive definitive operator 
—xD^ = aiid must not be confused with a power of the first order 

differential operator = for which the symbol is —in . An alternative 
illuminating notation for the symmetric fractional derivative is due to 
Zaslavsky|50]. and reads 

(3.5) 

For < a < 2 and |^^| < min {a, 2 — a} the Riesz-Feller derivative reads 
.D, m ^ ^ {.n lia . .)./2] f ^^^^^^ 



M ■ w m /9i r fix -0- fix) ,A 



(3.6) 



3.2 The Caputo fractional derivative 

For a sufficiently well-behaved function f{t) we define the Caputo time- 
fractional derivative of order /3 with < /3 < 1 through its Laplace transform 

£ {,D^/(t); 4 =//(.) -/^V(0+), 0</3<l. (3.7) 
This leads us to define 

./^f/(t):= r(W)X7^^' (3.8) 
i|/(t), ^^1. 

For the essential properties of the Caputo derivative, seejH [131 SH]- 
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3.3 The fundamental solution of the space-time frac- 
tional diffusion equation 

Let us note that the solution u{x, t) of the Cauchy problem (3.1)-(3.2), known 
as the Green function or fundamental solution of the space-time fractional 
diffusion equation, is a probabihty density in the spatial variable x, evolving 
in time t. In the case a = 2 and /3 = 1 we recover the standard diffusion 
equation for which the fundamental solution is the Gaussian density with 
variance cr^ = 2t. Sometimes, to point out the parameters, we may denote 
the fundamental solution as 

u{x,t) = G%{x,t) , (3.9) 

For our purposes let us here confine ourselves to recall the representation in 
the Laplace-Fourier domain of the (fundamental) solution as it results from 
the application of the transforms of Laplace and Fourier to Eq. (3.1). Using 
S{k) = 1 we have: 



hence 

■u(k, s) = ^(k, s) = — 5 . (3.10) 

^ ' ^ a./iv ' ^ _^ l^ja^fc'signK ^ ^ 

For explicit expressions and plots of the fundamental solution of (3.1) in the 
space-time domain we refer the reader to Mainardi, Luchko and Pagnini[29]. 
There, starting from the fact that the Fourier transform of the fundamental 
solution can be written as a Mittag-Leffler function with complex argument. 



u 



:K,t)=Gl,{K,t) = E,[-\K\^^^^^S^^t^) , (3.11) 



these authors have derived a Mellin-Barnes integral representation of 
u{x,t) = G^ij{x,t) with which they have proved the non-negativity of the 
solution for values of the parameters {a, 6, /5} in the range (3.2) and analyzed 
the evolution in time of its moments. The representation of u{x,t) in terms 
of Fox if-functions can be found in Mainardi, Pagnini and Saxena[3T]. see 
also Chapter 6 in the recent book by Mathai, Saxena and Haubold|34]. 
We note, however, that the solution of the STFD Equation (3.1) and its 
variants has been investigated by several authors; let us only mention some 
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of themP ElElinilSlE^lESlESlESlESlliniilllMlll^ where the 

connection with the CTRW was also pointed out. 

In particular the fundamental solution for the space fractional diffusion {0 < 
a < 2, /3 = 1} is expressed in terms of a stable density of order a and 
skewness 9, 

=ri/°L^(x/t^/°), -oo<x<+cx), t>0. (3.12) 

whereas for the time fractional diffusion {a = 2, < /3 < 1} in terms of a 
(symmetric) M- Wright function of order /3/2, 

G%{x, t) = ^t-^/^ Mp/2 {\x\/t^^'^) , -oo < X < +00 , t > , (3.13) 

For the standard diffusion {a = 2 , /3 = 1} we recover the Gaussian density 

GUx,t) = -^e-^V(4t) = r 1/2^0^^/^1/2^ ^ It-^'M,/, ilxl/t^^) . 
Zynt / 

Let us finally recall that the M- Wright function does appear also in the 
fundamental solution of the rightward time fractional drift equation, 

d 

tD^u{x,t) = -—u{x,t) , -cx)< X < +CX) , t > . (3.14) 
Denoting by Q^{x,t) this fundamental solution, we have 

g;ix,t) = h~'''^{¥)^ (3.15) 

[0, a; < 0, 

that for /3 = 1 reduces to the right running pulse 6{x — t) for a; > 0. For 
details see [151 150] . 

3.4 Alternative forms of the space-time fractional 
diffusion equation 

We note that in the literature there exist other forms alternative and 
equivalent to Eq. (3.1) with initial condition u{x,0) = uo{x) including the 
case Uo{x) = 6{x). For this purpose we must briefly recall the definitions of 
fractional integral and fractional derivative according to Riemann-Liouville. 
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The Riemann-Liouvillc fractional integral for a sufficiently well befiaved 
function f{t) {t > 0) is defined for any order // > as 



{t-Tr-'f{r)dT. 



(3.16) 



m Jo 

We note the convention = I (Identity) and the semigroup property 

,J'^,r^,J\J''^,J'^+'', ^^>0,u>0. (3.17) 

The fractional derivative of order n > in the Riemann-Liouville sense is 
defined as the operator tD'^ which is the left inverse of the Riemann-Liouville 
integral of order // (in analogy with the ordinary derivative) , 



fi>0. 



(3.18) 



If m denote the positive integer such that such that m — 1 < /j, < m , we 
recognize from Eqs. (3.16) - (3.18): 



tD^^fit) := tD^tJ^-^f{t) 



(3.19) 



Then, restricting our attention to a order /3 with < /3 < 1 (namely m — 1) 
the corresponding Riemann-Liouville fractional derivative turns out 



dt 
d_ 

dt 



' f{r)dT 
r(l-/3)yo {t-ry. 



< /3 < 1, 

^ = 1. 



(3.20) 



Then we get the relationship among the Caputo fractional derivative with 
the classical Riemann-Liouville fractional integral and derivative: 



tD^.m := tJ'-^D^m = [fit) - /(O)] = ,D^fit) 



/(O) 



(3.21) 

and, as a consequence, the equivalence of (3.1) with the following problems 

u{x, t) - u{x, 0) = tJ'^ :cDq u{x, t) , u{x, 0) = 5{x) , (3.22) 
d 



dt 



u{x, t) = tD^-^ :cDe u{x, t) , u{x, 0) = 5{x) . 



(3.23) 
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4 Analytic and stochastic pathways to sub- 
ordination in space-time fractional diffu- 
sion 



Our starting key-point to introduce the analytical and stochastic approaches 
to subordination in space-time fractional diffusion processes is the fundamen- 
tal solution of the space-time fractional diffusion equation in the Laplace- 
Fourier domain given by (3.10). 



4.1 The analytical interpretation via operational time. 

Separating variables in (3.10) and using the trick to write l/{z + a) for Re(2;-|- 
a) > as a Laplace integral 



z + a 



dp 







we have, identifying p := as operational time, the following instructive 
expression for (3.10): 

^{k,s) = J^"[exp (-t*\K\"i^^'^'^'^)] [s^-^exp {-t,s^)] dt, . (4.1) 

We note that the first factor in (4.1) 

:= exp (-i,|Ac| V«ig^«^) (4.2) 

is the Fourier transform of a skewed stable density in x, evolving in 
operational time i*, of a process x — y{t*) along the real axis x happening 
in operational time i*, that we write as 

U,eM ^t:'/'^Li{x/tl/'^) . (4.3) 

We can interpret the second factor 

g^(t,, s) := s'^-^exp i-ty) (4.4) 

as Laplace representation of the probability density in evolving in i of a 
process = t*(t), generating the operational time from the physical time 
t, that is expressed via a fractional integral of a skewed Levy density as 

q0{t*,t) = t^'/f' tJ'-^ L-/{t/t\/P) = t-^ M^itjt^) , (4.5) 
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see Eq. (2.26). To prove that qp{t^,t) (surely positive for > 0) is 
indeed a probability density we must further prove that is normalized, 

roo 

/ C[i3{t*, t) dt^ — 1 . For this purpose it is sufficient to prove that its Laplace 

transform with respect to is equal to 1 for = 0. To get this Laplace 
transform qp{s^,t) wc proceed as follows. Starting from the known Laplace 
transform with respect to 

qp{U, s) = s^-^exp i-UsI^) , (4.6) 

we apply a second Laplace transformation with respect to with parameter 
to get 

so, by inversion with respect to t 

POO 

M^*:t)= / e-'*'q0{U,t)dU = Ef){-sj'^), (4.8) 
and setting s* = 

poo 

q^iU,t)dU^E0{O)^l. (4.9) 



i. 



Weighting the density of x = y{t*) with the density of t^, = t^{t) over < 
t < oo yields the density u{x,t) in x evolving with time t. 

In physical variables {x,t}, using Eqs. (4.1) ~ (4.5), we have the 
subordination integral formula 



u[x 



POO 

fa,e{.x,t^)qf}{t^,t)dt^, (4.10) 
it,=o 

where fa,e{x,t^) (density in x evolving in t^) refers to the process x = y{t^) 
(i* — > x) generating in "operational time" the spatial position x, and 
qp{t*,t) (density in t^ evolving in t) refers to the process = t*(i) {t — >■ t*) 
generating from physical time t the "operational time" i*. 

Our aim is to construct a process x = x{t) whose probability density is u{x^ t), 
density in x, evolving in t.. We will soon find justification for denoting the 
variable of integration by t*. We will exhibit it as the "operational time" for 
our fractional diffusion process, and for distinction we will call the variable 
t its "physical time" In fact fa,eix,t^) is a probability density in x e IR , 
evolving in operational time > and qp{t^,t) is a probability density in 
t^ > 0, evolving in physical time t > 0. 
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4.2 Stochastic interpretation. 

Clearly fa,e{x,t^) characterizes a stochastic process describing a trajectory 
X — y{t^:) in the {t^,x) plane, that can be visualized as a particle travelling 
along space x, as operational time is proceeding. Is there also a process 
t* = t*(t), a particle moving along the positive t^, axis, happening in physical 
time t? Naturally we want t^,{t) increasing, at least in the weak sense, 

^2 > ^1 =^ t*{t2) > t*{ti) . 

We answer this question in the affirmative by showing that, by inverting 
the stable process t = whose probability density (in t, evolving in 

operational time t*) is the extremely positively skewed stable density 

rp{t,Q^t:'/^L/{t/tl/^). (4.11) 

In fact, recalling 

7p{s,Q^eicp{-ty), (4.12) 

there exists the stable process t = t{t^,), weakly increasing, with density in t 
evolving in given by (4.11). We call this process the leading process. 

Happily, we can invert this process. Inversion of a weakly increasing 
trajectory means that horizontal segments are converted to vertical segments 
and vicevcrsa jumps (as vertical segments) to horizontal segments (in 
graphical visualization). 

Consider a fixed sample trajectory t = t(t*) and its also fixed inversion 
t* = t*(t). Fix an instant T of physical time and an instant of operational 
time . Then, because t — t{t^) is increasing, we have the equivalence 

t,{T) <T,^T< t{T,) , 

which, with notation slightly changed by 

U{T) ^t[, n^U, T ^t, t{n) t' , 

implies 

/ g«,tX= / rp{t\U)dt\ (4.13) 
Jo Jt 

for the probability density q{t*,t) in evolving in t. It follows 

Q rOO POO Q 

q{t*, ^) = gfj^ t*) dt' = gf^0it'^ Q dt' ■ 
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We continue in the s*-Laplace domain assuming t > 0, 



g(S*,t) = {s^rp{t\s^) - 5{t')) dt' . 

It suffices to consider t > 0, so that we have 5{t') = in this integral 
Observing from (4.12) 

= ^, (4.14) 



s* + 



we find 
so that 

finally 



7p{t,s,) = l3t^-^E'p{-s.t''), (4.15) 

/oo 
s,(5t'^~^ E'p{-sJ^) dt' = Efsi-sJ^) , (4.16) 



q{t.,t)=t-^Mp{t,/t^), (4.17) 
From (4.16) we also see that 

?(s*,s) = ^ = g^(s*,s), (4.18) 

+ 

implying (4.6) and, see (4.7), 

q{t,,t) = qp{t.,t) , (4.19) 

so that indeed the process = t*(t) is the inverse to the stable process 
t = t(t^:) and has density qp{t^,t). 

Remark. The process at hand, t* = t*(t), which is referred to as the inverse 
stable subordinator, is honoured with the name " Mittag-Leffier process" by 
Meerschaert at al. EZ]- Honouring this process by the name of Mittag- 
Leffier can be justified by the fact that by (4.16) the Laplace transform 
of its density is a Mittag-Leffier type function or by the fact that it is a 
properly scaled diffusion limit of the counting function N{t) of the fractional 
generalization of the Poisson process whose residual waiting time probability 
is the Mittag-Leffier type function Epi^—t^), see recent papers of ours[TT| [T3]. 
In view of its probability density it may also be called the M-Wright process. 

Stipulating that there exists a weakly increasing process = t*(t) with 
density g^(t^.,t) we can analogously find the density of its inverse t = t(t^,) 
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which comes just as r^(t, However, in the context of our here presented 
considerations not being allowed to know that such process = t*(t) exists, 
we have taken as a gift from God the process t = t{t*) and shown by its 
inversion that there exists a process t* = with the desired properties. 

From the density rii{t,t^) of the leading process t = we have found the 
density of the directing process as given by the Laplace transform 

pair (4.6), that is 

g/}(t*, t) qj3{t^, s) = s^"^ exp (-t*s^) . 

In physical coordinates we have (4.5) and (4.17), so also an expression 
through an M- Wright function, 

qpiU, t) = tJ'-^ rp{t, t.) = Mp{t./t^) , (4.20) 

see Eq. (2.26). 



4.3 Evolution equations for the densities rp{t^t^) oi t = 

t{U) and qpit^.t) of = t^{t). 

The Laplace-Laplace representation of the density Tjjit^t^) of the process 
t — t{t^) is, according to (4.14), 

F,(.,..) = ^. 

This imphes 

s* r^{s, s*) - 1 = -s'^ r^{s, s*) , 

and by inverting the transforms and observing the initial condition r/3(t, = 
0) = 6{t) we arrive at the Cauchy problem 

J-r^(t, Q = - tD^.rpit, t.) , rp(t, = 0) = 5(t) . (4.21) 

Because it suffices to consider only f > where 5{t) — 0, we need not 
introduce a singular term on the right-hand side. 

The Laplace-Laplace representation of the density qp{t^,t) of the process 
U — t*{t) is, according to (4.18), 



(s*,s) 



+ 
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This implies 

s) — ^ — — s) , 

and by inverting the transforms and observing the initial condition {t^:,t — 
0) — we arrive at the Cauchy problem 

tD^^qpiU, t) = --^Qisit*, t) q^iU, t = 0) = 6{Q . (4.22) 

Because it suffices to consider only > where S{t^) — 0, we can ignore the 
delta function on the right-hand side. 

Remark The fractional differential equations in the above Cauchy problems 
have the same form. By replacing t by and r by g one of them goes over into 
the other. However, in the first problem the delta initial condition refers to 
the fractional derivative (of order f3), in the second problem to the ordinary 
(first order) derivative. These equations are akin with the time-fractional 
drift equation treated in (3.14) and (3.15), with different coordinates and 
proper initial conditions, as explained above. The process t = t{t^) of the ffist 
problem is a positive-oriented (extreme) stable process, whereas the process 
is a fractional drift process, see (3.14)-(3.15) with x replaced by 
t*. The reason for the two evolution equations to have the same form is 
that the described two processes are inverse to each other, their graphical 
representations coincide just by interchanging the coordinate axes. The delta 
initial condition for each equation is given at value zero of the evolution 
variable for the variable in which the solution is a density. 

4.4 The random walks. 

We can now construct the process x — x{t) for the position x of the particle 
depending on physical time t as follows in two ways. With the variable t 
(physical time), (operational time), x (position), we have the processes 
(i), (ii) and (iii), as follows: 

(i) t = t{t^) with density Vjiit^t^) in t, evolving in t^, the leading process, 

(ii) x = y{t^) with density fa,e{x,t^) in x, evolving in t^, the parent process, 

(iii) — t^{t) with density qp{t^,t) in t^, evolving in t, the directing process. 
Observing that the processes (i) and (iii) are inverse to each other, and 
taking account of the subordination integral (4.10), we define the space-time 
fractional diffusion process as the subordinated process 

X = x{t) = y{t,{t)) . (4.23) 
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Simulation of a trajectory for the subordinated process means: generate in 
running physical time t the operational time t^,, then the operational process 

y{t*)- 

Now, the Mittag-Leffler (or M- Wright) process t* = t*(t) is non-Markovian 
and not so easy to simulate. The alternative (we call it "parametric 
subordination") is to produce in dependence of the operational time the 
processes (i) and (ii) and then eliminate from the system 

t = t{Q , X = y{Q , (4.24) 

to get X = x{t) from x = y(t^) by change of time from t^. to t. 

We can produce a sequence of precise snapshots oit = t{t^.) and x = y(t^) in 
the (t*,t) plane and the (t*,x) plane by setting, with a step-size > 0, 

t*,n = , tn = T^^i + T*^2 ■ ■ ■ + T^^n , Xn = X^:^i + X^^2 + " • " + ^*,ra 5 (4-25) 

taking for k = l,2,...,n each T^.^^ as a random number with density 
T^^^^ L^^ (t/rl^^) and each X^^^k as a random number with density 
r* ^^'^ L^^{x/tI^°') , corresponding by self-similarity to the step t^. 

We can do this by taking random numbers Tk and Xk with density L^^(t) 

and L^{x), respectively, and then with r = r*^^ and h = r^", setting 
T^^^k = TTk, X^^k = hXk- In other words: we produce (a renewal process 
at equidistant times with reward) a positively oriented random walk on the 
half-line t > and a random walk on — oo < a; < +00 with jumps at 
equidistant operational time instants t^, = nr^,. We recognize the scaling 
relation jh'^ = 1, analogous to that used by us in earlier papers of ours on 
well-scaled passage to the diffusion limit in CTRW under power law regime, 
see [m [151 [16]. Methods for producing stable random deviates can be found 
in the books [25] 

Finally, we transfer into the (t,x) plane the points with coordinates t„,x„ = 
yn and so obtain a sequence of precise snapshots of a true process x = x(t). 
Finer details of the process x = x{t) become visible by using smaller values 
of the operational step-size r*. 
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5 Graphical representations and Conclusions 



We recall that, denoting the physical time with t, the operational time with 
t*, the physical space with x the density of the fractional diffusion process 
turns out to be given by the following subordination integral, see (4.10), 

POO 

u{x,t) = / fcfi{,x,U)qii{U,t)dU, (5.1) 

where fa^e{x,t^,), is the density (in x evolving in t^) of the parent process 
X = y{t^) = a;(t(t*)) and qpit^^.t) is the density (in evolving in t) of the 
directing process 

By using the Fourier-Laplace pathway we recall the two densities related to 
the parameters a, 9, (3 from (4.3) - (4.5), 

Ue{^,Q=t:'/^Li{x/t'/"), (5.2) 

q0{t*,t) = t-^lf" tJ''^ L/ [t/tl/^] = t-^ Mp{U/tP) , (5.3) 

where L refers to the Levy stable density and M to the Wright function, both 
introduced in Section 2. But for the parametric subordination the relevant 
density is ri^{t,t^,) governing the leading process, a density in the physical 
time t evolving with the operational time t^,: it turns out to be the unilateral 
Levy density of order (3, namely, see (4.11), 

r,it,t.)=t;'/n/{t/tl/^) . (5.4) 

We have shown in Section 4.4 that in our approach (referred to as parametric 
subordination) the process of space-time fractional diffusion (non-Markovian 
for P < 1) can be simulated by two Markovian processes governed by 
stable densities, provided by fa,e{,x,t^) and r^it^t^), as pointed out in our 
2007 paper with Vivoli[18j, where we have dealt with the CTRW model. 
There, before passing to the diffusion limit, we have two Markov processes 
happening on a discrete set of equidistant instants (for simplicity the non- 
negative integers) n = 0,1,2,..., meaning = 1, one of them moving 
randomly rightwards along t > 0, the other moving randomly on the real 
line —oo < x < oo with jumps T„ and Xn, respectively, for n > 1. By 
summing the "waiting times" and the jumps from 1 to we obtain 
sequences of jump instants t = tn and positions x = x„,that we display in the 
[t, x) plane. In fact, CTRW is the virgin form of parametric subordination. 
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We note that our approach is akin to that based on two stochastic differential 
equations, known in Physics as Langevin equations, see|71 ES]- We have 
indicated these two stochastic differential equations in|T8]. Here now we 
content ourselves with referring to the above cited papers. 

In Section 4.4 we have splitted the fractional diffusion process into three 
processes (i), (ii), (iii), each of them containing two of the three coordinates: 
space X, physical time t, operational time t*. We simulate the leading process 
by a random walk [rwi), the parent process by a random walk {rw2), and the 
subordinated process (which yields the desired trajectory) by a random walk 
{rw). The inversion of (rwi) gives us a random walk {rw^j for simulation of 
the directing process = t^,{t). Essentially, we need to carry out only (rwi) 
and {rw2) according to the equations (4.25). By transferring the points 
(t„, Xn) into the {t, x) plane we get the random walk (rw) as visualization 
of a random trajectory x = x(t) = y(t^(t)) according to the subordinated 
process which is our space-time fractional diffusion process of interest. 

To make transparent the situation we display as a diagram in Fig. [2] the 
connections between the four random walks. It is now instructive to show 




Figure 2: Diagram for the connections between the four random walks 
(rwi), {rw2), (rws) and (rw), related to the leading, parent, directing and 
subordinated processes, respectively. 

some numerical realizations of these random walks for two case studies of 
symmetric {9 = 0) fractional diffusion processes: {a = 2, (3 = 0.80}, {a = 
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1.5, /3 = 0.90}. As explained in a previous sub-section, for each case we 
need to construct the sample paths for three distinct processes, the leading 
process t = t(t*), the parent process x = y{t^) (both in the operational time) 
and, finally, the subordinated process x = x{t), corresponding to the required 
fractional diffusion process. 

We shall depict the above sample paths in Figs. |3l HI |5] respectively, devoting 
the left and the right plates to the different case studies. 

Plots in Figj3] (devoted to the leading process, the limit of (rwi)) thus 
represent sample paths in the (t*, t) plane of unilateral Levy motions of order 
p. By interchanging the coordinate axes we can consider Fig|3]as representing 
sample paths of the directing process, the limit of {rw^). 

Plots in FigJH (devoted to the parent process, the limit of {rw2)) represent 
sample paths in the (t*,x) plane, produced in the way explained above, for 
Levy motions of order a and skewness 6 = (symmetric stable distributions). 

By the indicated method, see (4.25), we have with (for simplicity) r^, = 1, 
^ = (symmetry) produced 10000 numbers tn and corresponding numbers 
yn- Plotting the points (t*,n,t„) into the (t^,t) (operational time, physical 
time) plane, the points (t^:^n,yn) into the (t*,x) (operational time, position 
in space) plane we get Figj3]and FigJUfor visualization of (rwi) and {rw2), 
respectively. Fig. [5] (as a visualization of (rw) is obtained by plotting the 
points (tn,yn) into the (t,x) (physical time and space) plane. 

Actually, we have invested a little bit more work in producing the figures. 
Namely, to make visible the jumps as vertical segments, we have in Fig|3] 
connected the points ((t*,n,^n) and ((t*,n+i,^n) by a horizontal segment, the 
points ((t*,n+i,^n) and {(t^^n+i,in+i) by a vertical segment. Analogously in 
FigJUwith the indexed i replaced by indexed y. In Figj5]we have connected 
the points {in, yn) and (tn+i, Vn) by a horizontal segment, the points (tn+i, yn) 
and (tn+i,yn+i) by a vertical segment. 

Resuming, we can consider Figj3] as a representation of {rwi) for the 
leading process, or by interchange of axes as one of [rws) for the directing 
process, FigJH as one of {rw2) for the parent process, and finally FigE] as a 
representation of (rw) for the subordinated process which is our space-time 
fractional diffusion process. 
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Figure 3: A sample path for {rwi), the leading process t = 
LEFT: {a = 2 , ^ = 0.80}, RIGHT: {« = 1.5 , ^ = 0.90}. 
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Figure 4: A sample path for {rw2), the parent process x = y{t^). 
LEFT: {a = 2, P = 0.80}, RIGHT: {a = 1.5 , /3 = 0.90}. 





Figure 5: A sample path for (rw), the subordinated process x — x(t). 
LEFT: {a = 2 , ^ = 0.80}, RIGHT: {« = 1.5 , ^ = 0.90}. 
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0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

Figure 6: A sample path for the leading process t = t{t,,). 
LEFT: {/3 - 0.9, TV = 10^}, RIGHT: {/3 = 0.8, N = 10^}. 
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Figure 9: A sample path for the parent process x = y{t^). 
LEFT: {a = 2, N = 10^}, RIGHT: {a = 1.5, TV = 10^}. 
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Figure 10: A sample path for the parent process x = y{tt.). 
LEFT: {a = 2, N = 10^}, RIGHT: {a = 1.5, N = 10^}. 





Figure 11: A sample path for the parent process x = y{t^). 
LEFT: {a = 2, iV = lO^}, RIGHT: {a = 1.5, N = 10^}. 
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Figure 12: A sample path for the subordinated process x = x{t). 
LEFT: {a = 2, /3 = 0.80, N = 10^}, RIGHT: {a = 1.5, /3 = 0.90, N = 10^}. 
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Figure 13: A sample path for the subordinated process x = x{t). 
LEFT: {a = 2, /3 = 0.80, N = 10^}, RIGHT: {a = 1.5, /3 = 0.90, N = 10^}. 
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Figure 14: A sample path for the subordinated process x = x{t). 
LEFT: {a = 2, /3 = 0.80, N = lO^}, RIGHT: {a = 1.5, /3 = 0.90, N = 10^}. 
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We conclude by including additional Figures, showing the effect of taking 
smaller step-sizes r^,, equivalently larger values N of steps following our 
analysis for the CTRW^. Figs J6]l71l8l Figs MTOllTTl Figs HaTSUl show the 
effect of making the operational step-length r* smaller or, equivalently, the 
number N of operational steps larger for the sample paths of the leading, 
parent and subordinated processes, respectively. In these pictures = 
and we have taken = 10, = 100 and = 1000. Finer details will 
become visible by choosing in the operational time the step-length r* 
smaller and smaller. In the graphs we can clearly see what happens for finer 
and finer discretization of the operational time t^,, by adopting 10^, 10^, 10^ 
of number of steps. As a matter of fact there is no visible difference in the 
transition for the successive decades 10*^, 10^, 10^ of number of steps as the 
great majority of spatial jumps and waiting times even for very small steps 
r* of the operational time. This property also explains the visible persistence 
of large jumps and waiting times. 
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